A MATHEMATICAL MODEL OF SAP EXUDATION IN MAPLE TREES 
GOVERNED BY ICE MELTING, GAS DISSOLUTION AND OSMOSIS* 
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Abstract. We develop a mathematical model for sap exudation in a maple tree that is based on a purely 
physical mechanism for internal pressure generation in trees in the leafless state. There has been a long-standing 
controversy in the tree physiology literature over precisely what mechanism drives sap exudation, and we aim to 
cast light on this issue. Our model is based on the work of Milburn and O'Malley [Can. J. Bot., 62(10): 2101-2106. 
1984] who hypothesized that elevated sap pressures derive from compressed gas that is trapped within certain wood 
cells and subsequently released when frozen sap thaws in the spring. We also incorporate the extension of Tyree 
[in Tree Sap, pp. 37-45, eds. M. Terazawa et al., Hokkaido Univ. Press, 1995] who argued that gas bubbles are 
prevented from dissolving because of osmotic pressure that derives from differences in sap sugar concentrations and 
the selective permeability of cell walls. We derive a system of differential-algebraic equations based on conservation 
principles that is used to test the validity of the Milburn-O'Malley hypothesis and also to determine the extent to 
which osmosis is required. This work represents the first attempt to derive a detailed mathematical model of sap 
exudation at the micro-scale. 
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1. Introduction. Sap flow in deciduous trees is driven during the growing season by the 
process of transpiration, in which evaporation of water in the leaves draws groundwater from the 
roots to the crown [T4J [28] . As much as 90% of water taken up by the tree is transpired into the 
atmosphere to generate the large pressure differential needed to overcome gravity; only 10% of the 
water is actually consumed in photosynthesis reactions in the leaves to produce the sugars needed 
for growth and other life processes. These sugars are transported by the sap in dissolved form 
back to the trunk and roots where they are either consumed immediately or else stored as starch 
for later use. Transpiration halts in late fall when the leaves drop and the tree enters its winter 
dormant phase, which lasts until the spring thaw. At that time, the stored starch reserves are 
released and provide the energy needed to initiate budding and leaf growth before photosynthesis 
and transpiration can begin again. 

In the sugar maple (Acer saccharum), a process known as sap exudation occurs in the time 
period between the dormant leafless phase and the active transpiration phase, when certain pro- 
cesses are triggered that build up positive sap pressure and convert stored starches into sugars. 
The stem pressure is large enough that when the tree is tapped, sap exudes naturally from the 
tap-hole. The mechanism that generates this internal pressure is unique to the sugar maple and 
a few other related species (black and red maple, birch, walnut, and butternut) whose sap sugar 
content is also significant although not typically as high as sugar maple. Raw maple sap typically 
contains 2 to 3% sugar by weight (primarily sucrose) and can be processed by repeated boiling to 
produce sweet syrup and other edible maple products. 

The scientific study of sap exudation has a long history that has seen several competing 
hypotheses proposed for the mechanism that drives sap flow in spring. These hypotheses can be 
roughly divided into two classes: "vitalistic" models that require some intervention by living cells; 
and "physical" models that rely on passive, physical effects. Early in the 20th century, Wiegand's 
experiments [32j showed that even though up to one-quarter of the total volume of sapwood or 
xylem is filled with gas, the expansion and contraction of this gas in response to temperature 
changes alone could still not account for the observed pressure changes. He concluded that some 
active process initiated by living cells must be responsible, a claim that was supported by the 
experiments of Johnson [7]. 

Later studies questioned this vitalistic hypothesis, for example Stevens and Eggert [23] who 
attributed the pressure generation in leafless maple trees to volume changes arising from freezing 
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and thawing of sap. This hypothesis was further explored in subsequent years [201 125] . culminating 
in the ground-breaking paper of Milburn and O'Malley [T2] who proposed that gas trapped in the 
xylem is compressed by growth of ice crystals and subsequent uptake of water when the tree freezes 
in fall; during the spring thaw, they argued that this compressed gas generates the pressures that 
drive exudation. This freeze/thaw hypothesis was later supported by experimental evidence [8lll3] 
showing that gas expansion due to thawing sap is the primary cause of exudation pressure; they 
also found no evidence to suggest that living cells play a significant role. 

However, the Milburn-O'Malley model was unable to account for all observations, and so some 
authors [22j [26j [28] cast doubt on the ability of gas expansion alone to explain exudation, arguing 
that heightened pressures in the xylem will dissolve any gas bubbles within a few hours. Several 
studies [SJ [B] suggested that osmotic pressure, generated by differences in sugar concentration 
across semi-permeable structures that separate various cells in the xylem, could play a significant 
supporting role in sap exudation. Most notably, Tyree proposed [26J that a combination of the 
freeze/thaw mechanism and osmosis is the key to maintaining stable gas bubbles in the xylem over 
observed exudation time scales. 

This issue remains controversial and the current situation is perhaps best summed up by 
Ameglio et al. who concluded that "no existing single model explains all of the xylem winter 
pressure data," so that both vitalistic and physical effects must play a role [2] . The lack of a clear 
consensus in the tree physiology literature regarding the precise mechanisms driving sap exudation 
is compounded by the fact that no detailed mathematical model has yet been developed for the 
freeze/thaw mechanism, not to mention the other effects of osmosis and gas dissolution. 

The aim of this paper is therefore to develop a mathematical model that captures the essential 
processes driving sap exudation. We begin in Section [2] by providing details of tree physiology, sap 
flow, and the Milburn-O'Malley and Tyree hypotheses for pressure generation, at the core of which 
is the special role played by two types of non-living xylem cells called vessels and fibers. In Section[3] 
we derive a system of differential equations governing the multiphase gas/liquid/ice system that 
incorporates the effects of thawing sap in the fibers, dissolving gas bubbles in the vessels, and 
the osmotic pressure pressure gradient between the two. Numerical simulations in Section [4] are 
used to validate the model and conclusions are drawn regarding the various mechanism for sap 
generation in Section [5] 

2. Basic tree physiology and the Milburn O'Malley hypothesis. We first provide 
some background information on tree structure and the hydraulics of sap flow that can be found in 
texts such as [TU |2S]. We also describe the freeze/thaw model proposed by Milburn and O'Malley, 
and indicate how other mechanisms such as gas dissolution and osmotic pressure enter the picture. 

The Milburn-O'Malley model is a purely physical one that operates at the cellular level and 
depends fundamentally on the structure of the non-living cells making up the xylem. As pictured 
in Figure POT a), the xylem consists of long, hollow, roughly cylindrical cells of two types: fibers, 
that are the main structural members in the wood; and vessels, that have a larger radius and 
form the primary conduit for transmitting sap. An individual vessel can extend for up to tens 
of centimeters vertically through the xylem, and is subdivided into elements that are connected 
end-to-end via perforated plates to form a continuous hydraulic connection. Vessels are connected 
to each other via numerous pits or cavities that perforate the vessel walls. The fibers are in turn 
subdivided into two sub-classes: libriform fibers and tracheitis. According to Cirelli et al. [5], 
the tracheid walls also contain pits through which they exchange sap with the vessels and hence 
integrating the tracheids into the water transport system of the tree. We will therefore ignore the 
tracheids, treating them as a part of the vessel network, and focus instead on the fibers. 

In comparison with tracheids, the libriform fibers (or simply "fibers") lack the pits that provide 
a direct hydraulic connection to the vessels and so they are often considered to play a primarily 
structural role. However, the recent experiments of Cirelli et al. [5] suggest that the fiber sec- 
ondary wall has a measurable permeability to water. Other experiments show that under normal 
conditions, maple trees (and other species that exude sap in spring) are characterized by gas-filled 
fibers and water-filled vessels [HI [26] . This is in direct contrast with most other tree species whose 
fibers are completely filled with water. Hence, it is reasonable to suppose that the presence of gas 
in xylem fibers may be connected with the ability of maple trees to generate exudation pressure 
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in the leafless state. 

The mechanism proposed by Milburn and O'Malley [12], and discussed in complete detail 
in [35], can be summarized as follows: 

• During late fall or early winter when the tree freezes, ice crystals form on the inner wall 
of the fibers and the growing ice layer compresses the gas trapped inside. 

• In early spring when temperatures rise above freezing, ice crystals in the fiber melt and 
the pressure of the trapped gas drives the melted sap through the fiber walls into the 
vessel, hence leading to the elevated pressures observed during sap exudation. Tyree [25] 
estimates that this effect is responsible for a pressure increase of roughly 30-60 kPa. 

• Exudation pressure is enhanced by the gravitational potential of sap that was drawn into 
the crown during the previous fall, and once thawed is then free to fall toward the roots. 

As mentioned before, this model can account for the initiation of exudation pressure but it cannot 
explain how pressures are sustained over periods longer than about 12 hours because gas bubbles 
will dissolve when pressurized. Tyree |26j explains how surface tension effects require that any gas 
bubble is at higher pressure than the surrounding fluid. As a result, Henry's law necessitates that 
the concentration of dissolved gas in the region adjacent to the bubble is proportional to interfacial 
pressure. Since the gas concentration is highest there, dissolved gas diffuses away from the bubble 
and eventually causes the bubble to disappear. 

One clue to identifying a mechanism that can sustain gas bubbles over longer time periods is 
provided by experiments that show maple trees exude sap only when sugar is present [8]. With 
this in mind, Tyree |26) proposed a modification of the Milburn-O'Malley mechanism in which 
osmotic pressures arising from differences in sucrose concentration could permit gas bubbles to 
remain stable. In particular, he suggested that the fiber secondary walls are selectively permeable, 
permitting water to pass through but not larger molecules like sucrose. Therefore, the ice layer 
that forms on the inner surface of the fibers in fall/ winter is composed of pure water, which leads 
to an osmotic potential difference between the fibers and the vessels. 



3. Derivation of the model equations. We now derive a compartment model describing 
the dynamics of gas, water and ice in a single vessel that is in contact with the surrounding fibers. 




(a) 3D view of the xylem, showing vessel el- 
ements surrounded by fibers and illustrating 
the difference in the cell diameters. The pits 
forming the hydraulic connections between the 
cells are also shown. Source: J29f . 



(b) 2D cross-section highlighting the vessel elements ( ve), 
libriform fibers (l-f), and fiber tracheids (f-t). Source: 
Cirelli et al. Fig. 6b] (reproduced with permission, Ox- 
ford Univ. Press). 



Fig. 2.1. Two views of the layout of vessels and fibers in the xylem. 
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FlG. 2.2. Picture of ice and compressed gas taken from Milburn and O'Malley ]12\ Fig. 7] (© 2008 Canadian 
Science Publishing, reproduced with permission). We are concerned here only with the "warming sequence" in the 
bottom row, numbered 5-8. The fiber is the rectangular structure on the left of each image, while the vessel is 
represented by the vertical channel on the right (not drawn to scale). 



To begin with, we restrict ourselves to the time period before the ice in the fiber melts completely. 
The geometry is depicted in Figure IQ1 which is an idealized view of Figure l2.1f b). Before deriving 
the governing equations, we make the following simplifying assumptions: 

Al. We consider only the thawing phase, assuming that all water within the fibers is initially 
frozen. 

A2. Gravitational potential differences due to height can be neglected since we focus only 
locally on a small section of the xylem. 

A3. We consider a single vessel element and ignore interactions between adjacent elements. 
This is a reasonable assumption if all elements experience similar conditions, and any gas 
contained within an element is isolated by the perforated plates from that of its neighbors. 

A4. The vessel is surrounded by N identical fibers. Although our model considers only a single 
fiber explicitly, the fluxes and other fiber quantities are scaled by a factor N to determine 
the total influence of all N fibers on the vessel. 

A5. The fiber and vessel are hollow and cylindrical in shape, with interior radii R* and R v , 
and heights and V respectively. Both have large aspect ratio so that R/L <C 1. 

A6. We assume the problem is one-dimensional, with the domain extending horizontally from 
the center of the fiber to the rightmost outer edge of the vessel (pictured at the top of 
Figure 13.11) . The coordinate system is chosen so that the horizontal (x) axis is directed 
along the line joining the centers of the fiber and vessel and with origin x = at the center 
of the fiber. Referring to Figure [3TTI the fiber/vessel interface is located at x = R* and 
the center of the vessel at x = R? + R v . 

A7. Gas is present in both fiber and vessel. Although the existence of gas in the fiber is a 
fundamental assumption of the Milburn-O'Malley model, the gas in the vessel is a novel 
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feature introduced in our model. Because water is incompressible and the xylem is a 
closed system when the tree is frozen, pressure cannot be transferred between the fiber 
and vessel compartments unless gas is also present in the vessel. Further justification 
for this assumption is provided in the papers |22l 127] . which indicate that maple trees 
experience winter embolism (bubble formation) when gas in the vessels is forced out of 
solution upon freezing. 

A8. Gas in both fiber and vessel takes the form of a cylindrical bubble located at the center 
of the corresponding cell. This seems reasonable in the vessel where the surface tension 
is of the order of a/R v ~ 4 x 10 3 Pa, which is several orders of magnitude smaller than 
the typical gas and liquid pressures. In the fiber, the smaller radius gives rise to a much 
larger surface tension (a/Rf « 2 x 10 4 Pa) which though still small relative to gas and 
liquid pressures could still potentially initiate a break-up into smaller bubbles owing to 
the Plateau-Rayleigh instability [?]. However, regardless of the precise configuration of 
the gas in the fiber, we assume that the net effect on gas pressure is equivalent to that of 
a single large bubble. 

A9. Heat from outside the tree enters from the right in Figure l3~T1 Consequently, the sap in 
the vessel is taken to be initially in liquid form and the ice in the fiber begins melting on 
the inner surface of the fiber wall. 
A10. Gas and ice temperatures in the fiber can be taken as constant and equal to the freezing 
point. This is justified by the fiber length scale being so much smaller than that of the 
vessel (Rf < R v ). 

All. Time scales for heat and gas diffusion are much shorter than those corresponding to ice 
melting and subsequent phase interface motion (which our simulations show are on the 
order of minutes to hours). This is justified by considering the various diffusion coefficients 
(D) and a typical length scale (x w 10~ 5 m), and in each case calculating the corresponding 
diffusion time scale t ~ x 2 /AD: 

- diffusion in gas bubbles: using the air self-diffusion coefficient D s» 2 x 1CP 5 m 2 /s, 
the time scale is t ps 1.3 x 10 -6 s; 

- diffusion of dissolved gas: diffusivity D m 1 x 10~ 9 of air in water gives t ss 0.025 s; 

- diffusion of heat: thermal diffusivities D ss 1.9 x 10~ 5 (for air) and 1.4 x 10~ 7 (for 
water) yield t ~ 1.3 x 10~ 6 and 1.8 x 10~ 4 respectively. 

Therefore, heat transport is essentially quasi-steady and the temperature equilibrates 
rapidly to any change in local conditions. Also, gas concentration and density (in both 
gaseous and dissolved forms) can be taken as constants in space. 
A12. A significant osmotic pressure arises owing to differences in sugar content between liquid 
in the vessel and fiber. This is motivated by recent experiments which show that the fiber 
cell wall is permeable to water but not to larger molecules such as sucrose [5], suggesting 
that ice contained in the fiber is composed of pure water. A simple order of magnitude 
estimate based on a sap sugar concentration of 2% then implies that the osmotic pressure 
n ~ O(10 5 ) Pa, which is of the same order as typical gas and liquid pressures. 
A13. There is no need to consider modelling the system beyond the time when either fiber or 
vessel bubble is completely dissolved, since then no further exchange of pressure is possible. 
Model equations are derived for the fiber and vessel compartments in the next two sections. 
For easy reference, all parameters and variables are listed in Table [3~T1 along with physical values 
where appropriate. There are some similarities that can be found with other models developed 
for related problems such as ice lensing in porous soils |24j and freeze/thaw processes in the food 
industry [11] : however, neither problem contains all of the physical mechanisms encountered during 
sap exudation. 

3.1. Fiber model. We now develop a set of governing equations for the fiber that is restricted 
to the time period when ice is present, so that an ice layer separates the gas from the water adjacent 
to the wall. The modifications that are required to the model once the ice layer has disappeared 
and the fiber bubble is in direct contact with the sap are described in Section 13.51 

The inner fiber wall is initially coated by an ice layer of uniform thickness, inside of which 
is a cylindrical volume of compressed gas with radius s g i(t). As the ice layer melts, an ice/water 
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Fig. 3.1. (Bottom) A 2D cross-section through a fibei — vessel pair showing the water, ice and gas regions, and 
the moving interfaces. (Top) The ID region corresponding to our simplified model geometry is outlined with a red, 
dashed box in the bottom figure. The locations of the various phase boundaries are indicated, with the origin x = 
at the center of the fiber. 



interface x — Si W (t) appears at the inner edge of the fiber wall (si w (0) — Rf) and propagates 
inward. The speed of the interface depends not only on the volume difference between ice and 
water but also on the seepage of liquid through the porous fiber wall. The gas pressure that drives 
water through the wall into the vessel is transmitted through the ice layer to the adjacent water. 

If we denote by Vj(t), Vjj(t) and V/(t) the volumes of gas, water and ice within the fiber, 
then owing to the cylindrical symmetry 

V g f (t)=irLf Sgi (t) 2 , (3.1) 

Vjj(t)=TrLf((Rf) 2 -s iw (tf), (3.2) 

V/(t) = xlS (s lw (t) 2 ~ s g i(t) 2 ) , (3.3) 

which also satisfy 

V g f (t) + Vf{t) + V/(t) = nLf (Rff =: V* . (3.4) 

Recalling Assumption lAlOl which takes the gas temperature Tf to be constant and equal to the 
freezing point T c — 273.15 °K, the gas density varies only due to changes in volume so that 



V f (0) 

The corresponding gas pressure is given by the ideal gas law 

fU\ P 9 ( 



_ p f g (t)nT c 



1 y 
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Table 3.1 

A list of nomenclature for all physical parameters and variables, including values in SI units. References are 
given for each numerical value, either as a citation to the literature (in square brackets) or else an equation number 
(in parentheses) for those values calculated from other parameters. 



Description | Units | Value | Ref. 



Constant parameters: 



A 


Area of vessel wall that is permeable 


m^ 


6.28 x 10" s 


(3.18J 


Cs 


Sucrose concentration in vessel sap 


mol/m 3 


58.4 






9 


Gravitational acceleration 


m/s 2 


9.81 




30] 


h 


Heat transfer coefficient 


W/m 2 °K 


10.0 




13 


H 


Henry's constant 


- 


0.0274 




M 


kg 


Thermal conductivity of air 


W/m °K 


0.0243 




UJ 


k 


Thermal conductivity of water 


W/m °K 


0.580 




Uj 


K 


Cell wall hydraulic conductivity 


m/s 


1.98 x 10~ 14 




M 


Lf 


Length of fiber 


m 


1.0 x 10" 3 






L v 


Length of vessel element 


m 


0.5 x 10" 3 


: 


m 


Mg 


Molar mass of air 


kg/mol 


0.0290 




30] 


Ms 


Molar mass of sucrose 


kg/mol 


0.3423 




30] 


N 


Number of fibers per vessel 




16 




V 


Pressure scale 


N/m 2 


1.01 x 10 5 


|3.29| 


Rt 


Interior radius of fiber 


m 


3.5 x 10" 6 




33] 


R v 


Interior radius of vessel 


m 


2.0 x 10" 5 




33] 


K 


Universal gas constant 


J/mol °K 


8.314 




30] 


i 


Time scale 


s 


46.1 




T c 


Freezing temperature of ice 


°K 


273.15 






T a 


Ambient temperature of gas 


°K 


0.005 + T c 






V f 


Volume of fiber 


m 3 


3.85 x 10~ 14 


J3741 


yv 


Volume of vessel 


m 3 


6.28 x 10~ 13 


)3.28( 


w 


Thickness of fiber + vessel wall 


m 


3.64 x 10~ 6 


m 


A 


Latent heat of fusion 


J/kg 


3.34 x 10 5 




m 


pi 


Ice density 


kg/m 3 


917 






Pw 


Water density 


kg/m 3 


1000 




M 


p 


Density scale 


kg/m 3 


1.29 




26] 




Air-water surface tension 


N/m 


0.0756 




30 



Dependent and independent variables: 



c 9 


Dissolved gas concentration 


mol/m 13 


m 


Mass 


m 


V 


Pressure 


N/m 2 


r 


Radius of gas bubble in vessel 


m 


s 


Fiber interface location 


m 


t 


Time 


s 


T 


Temperature 


°K 


U 


Volume of melted ice 


m 3 


V 


Volume 


m 3 


X 


Spatial location 


m 


P 


Density 


kg/m 3 


Subscript (denotes phase): 




9 


Gas phase 




i 


Ice phase 




w 


Water phase 




Superscript (denotes cell component): 




f 


Fiber 




V 


Vessel 





where M g is the molar mass of air. This can be combined with (|3.1I) and (|3.5I) to obtain 



P f S)=V f M (?#) • (3-6) 



(*) 



As ice thaws, some volume of melt-water (denoted U{t)) is forced through the fiber wall into 
the vessel, while the remainder stays in the fiber sandwiched between the ice layer and the wall. 



8 M. CESERI AND J. M. STOCKIE 

Let the (constant) mass of gas in the fiber be m^ and the mass of water and ice be 

m{(t)= Pl v/(t), (3.7) 

ml{t) = p w v£(t), (3.8) 

where p w and pi are the densities of water and ice respectively. Conservation of mass in the fiber 
implies that 

m f g0 + m{(t)+m f w (t)+p w U(t)=m f , (3.9) 
where initially U(0) — and ml = m,g + mj(0) + to£,(0). Differentiating (|3.9p yields 

m{(t)+m{ u (t)+p w U(t) = 0, 



where the "dot" denotes a time derivative. Substituting (I3.4[) and (|3.7[) (|3.8[) into the last expres- 
sion yields 

- Pi V g f + (p w - Pi )vf + pjj = 0, 
and the volume terms can be replaced using (|3.1[) and (|3.2j) to obtain 

2 P iSgiS gl + 2(p w - pi)s iw s iw - -^jj [7 = 0. (3.10) 

Recall that the temperature in both gas and ice is assumed constant and equal to the freezing 
point T c . The water layer, on the other hand, is heated from the vessel side so that the water 
temperature T£(x, t) is not constant but instead obeys the steady-state heat equation and boundary 
condition 

d xx Tl = for all x € {s lw ,R f ) , (3.11a) 
Tl = T c atx = s iw . (3.11b) 

The speed Si W of the ice/water interface is determined by the Stefan condition 

^Pw^iw — kw^xTyj at X — Si w , (3.12) 

where A is the latent heat of fusion and k w is the thermal conductivity of water. 

In summary, (I3.6[) . (13.101) . (|3.11a[) and p. 121) represent a coupled system of equations for the 
fiber unknowns Si W (t), s g i(t), p-jjit), and T£(x,t). In the next section, we will derive an equation 
for U and also determine the missing boundary condition for T£ as a matching condition with the 
vessel temperature at the fiber/vessel wall x — . 

3.2. Vessel model and matching conditions. Let V^(t) and Vg{t) denote the volume of 
the water and gas compartments. The time derivative of the total volume (gas + water) must be 
zero so that 

±(v:(0) + NU(t) + V g v (t))=0, 

or simply 

Vg" = -NU. (3.13) 

This equation embodies the assumption that there is no flow between vessel elements so that the 
portion of the volume consisting of water equals the initial value V" (0) plus whatever water enters 
from the N surrounding fibers. Using the fact that the bubble radius and volume are connected 

by 

V?{i) = -KL°r{tf, (3.14) 
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equation (|3.13[) can then be rewritten as 



NU 

rf = -—^-. (3.15) 
2ttL v y ' 



Assuming that fibers are packed tightly around the vessel and that both vessel and fiber have wall 
thickness W , a simple geometric argument provides an estimate for the constant N: 

2R1 + W 1 ' 

We next focus on porous transport through the fiber wall which is governed by Darcy's law 

K A 

NU{t) = -— w (pl{t)-pf g {t)), (3.17) 

where p^(t) is the vessel water pressure, K is the hydraulic conductivity, W is the wall thickness, 
and A is the area of the wall that is permeable to water. In the absence of experimental measure- 
ments for hydraulic conductivity in maple, we use the value K = 1.98 x 10~ 14 m/s obtained for 
birch trees [16] . Finally, A is taken equal to the internal surface area of the cylindrical vessel 

A = 2irR v L v . (3.18) 

We next describe the temperature evolution within the various vessel compartments, recalling 
that the diffusive time scale is short in relation to either melting or gas dissolution. We therefore 
employ a quasi-steady approximation in which temperatures obey a steady-state heat equation. 
Denote by T g (x,t) the temperature within the gas bubble, and T^ 1 (a:,t) and T^ 2 (x, t) the water 
temperatures to the left and right of the bubble respectively. Then, in the water region on the left 
(adjacent to the vessel wall) T^^x, t) obeys 

dxxTZi = for a11 x e ( Rf , R f + R V -r). (3.19a) 

Both temperature and heat flux are continuous at the vessel wall 

T^=Tl a,tx = R f , (3.19b) 

k w d x T^ = k w d x Tl at x = R f , (3.19c) 

and at the water/gas interface 

T l wl = T v g at x = R f + R v - r, (3.19d) 

k w d x T^ = k g d x T v g at x = R f + R v - r. (3.19e) 

The gas temperature T g {x,t) inside the bubble obeys 

d xx T v g = for all x € (R f + R v - r, R f + R v + r) , (3.20a) 

with temperature and flux continuity conditions 

T; = T^ 2 at x = R f + R v + r, (3.20b) 

k g d x T v g = k w d x Tl 2 at x = R f + R v + r. (3.20c) 

Finally, the water temperature T^ 2 {x,t) to the right of the gas obeys 

d xx Tl 2 = for all x G (R f + R v + r, R f + 2R V ) , (3.21a) 
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with a Robin condition on the right-most boundary 

k w d x T^ 2 = h(T a - Tl 2 ) at x = R f + 2R\ (3.21b) 

where h is a convective heat transfer coefficient. The ambient temperature T a refers to the tem- 
perature in the neighboring layer of xylem cells; using a temperature difference of 10°K through 
the sapwood over a tap depth of 0.05 m, we obtain a rough estimate of T a ss 0.005°K for the 
temperature difference on the cell scale. 

The final component of the model is a description of the process whereby elevated pressures 
in the gas bubble cause the gas to dissolve at the gas/water interface. Because of the small bubble 
size, it is essential to take into account the effect of surface tension. Gas and liquid pressures in 
the vessel are connected by the Young-Laplace equation^, 

P v w (t) = p v g (x,t) - ^ atx = i?' + R V -r, (3.22) 

where a is the air-water surface tension and r(t) is the bubble radius. In this formula, the gas 
pressure can be written using the ideal gas law 

p v a (t)TZT^(x,t) 

M > ( 3 - 23 ) 

where p g (t) is the gas density which we have taken uniform in space according to Assumption I Al 11 
The concentration c v g {t) of dissolved gas in the sap is then related to the gas density by Henry's 
law 

c v g (t) = §-P v g W, (3-24) 

where H is a dimensionless constant. The gas density in the preceding equations decreases over 
time owing to dissolution and can be written 

pg(o)v fl "(o) - m 9 J g(g dv p g(o)v^(o) - M g cg(t)(^ - TO) 

where V g (t) is given by (13. 14|) and the integral in the first expression is taken over tl w , the annular 
portion of the vessel filled with liquid. 

In summary, equations (|3.15|) . (|3.17[) and (|3.19p (|3.25p represent a nonlinear system of differential- 
algebraic equations for the nine unknown functions r(t), p^(i), pZ(x, t), T%, x (x, i), T£, 2 (x, t), T g (x, i), 
p v Q (t), c " Q (t) and U (t) . The vessel and fiber solutions are coupled via matching conditions (13. 1T[) , 
([frtb]) and (j3~T9c)) . 

3.3. Osmotic effects. Tyree and co-workers have argued that a pressurized gas bubble in the 
xylem will dissolve entirely over a period of 12 hours or less, so that some other pressure-generating 
mechanism must operate to sustain the bubbles actually observed in xylem sap [5, 26, 28]. Recent 
experiments by Cirelli et al. [5] suggest that the fiber secondary wall is selectively-permeable, 
allowing water to pass through but preventing the passage of larger molecules such as sucrose 
and thereby generating a significant osmotic potential difference between fibers and vessels. This 
leads naturally to the hypothesis that osmotic pressure might provide the extra mechanism needed 
to enhance flow of sap from fiber to vessel and hence prevent bubbles in the fiber from totally 
dissolving once the ice layer has completely melted. 

Osmosis can be described mathematically by means of the Morse equation that relates osmotic 
pressure II within a solution to the dissolved solute concentration c s via II = c s TZT . Following the 



lr The Young-Laplace equation states that the pressure difference is proportional to the mean curvature, p g —p w = 
a y-f^ + wricrc the proportionality constant is the interfacial surface tension a. For a cylindrical gas bubble 

with radius r, the principal radii of curvature are Ri = r and R2 = 00. 
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approach in |18j . we apply this equation to the sap solutions in fiber and vessel and hence obtain 
a modified version of the melt volume equation (|3.17p 

K A 

NU = -T-nW ^ -P f g )-KTl(Rf,t) {cl - c{)] , 

where c{ and c v s are the sucrose concentrations in the fiber and vessel respectively. According to 
Tyree and Zimmermann [28], sucrose is present only in the vessels, so that c{ = and 

K A 

NU= -^W [(j>l-P S g )-KTl{Rf,t)cl\. (3.26) 
Pw9 

This equation replaces (|3.17|1 . but otherwise the governing equations remain unchanged. For the 
sake of simplicity, we assume that the vessel sucrose concentration c v s is a constant which for a 2% 
sucrose solution corresponds to taking c" = 0.02p w /M s s» 58.4 mol/m 3 , where the molar mass of 
sucrose is M s = 0.3423 kg/mol. This assumption can be justified by arguing that ray parenchyma 
cells (which provide the bulk of the sucrose to the tree vascular system) are so numerous in maples 
that they should be capable of maintaining the sucrose concentration at a relatively constant level. 

3.4. Non-dimensional variables. To simplify the governing equations, we reduce the num- 
ber of parameters in the problem by introducing the following change of variables: 

x = R v x, s = R v s, r = R v f, 

t = it, U = V V U, P = PP, (3.27) 

p = pp, c^j^-c, T — T c + (T a - T C )T, 

where an overbar denotes a dimensionless quantity. The melt volume is non-dimensionalized using 

V v =ttL v (R v ) 2 , (3.28) 
density and pressure scales are chosen based on initial values as 

p = pl(0) and p=pl(0) = ^, (3.29) 

and the time scale t will be specified shortly. 

The above expressions are then substituted into the dimensional equations from Section |3l 
after which we immediately simplify notation by dropping asterisks. The differential equations 
and boundary conditions p. lip and p.l9[) - (|3.21[) governing the temperatures T£(x,t), T^ ]1 (x,t), 
T£ 2 (x,t) and f%(x,t) become: 

Fiber water: dxxT£ = for all x £ (si W , S) , (3.30a) 

(3.30b) 

Vessel water (left): d^f^ = for all x € (6, 5 + 1 - f) , (3.31a) 

(3.31b) 
(3.31c) 
(3.31d) 
(3.31e) 

Vessel gas bubble: d xx T" = for all x G {5 + 1 - f, S + 1 + f) , (3.32a) 

(3.32b) 
(3.32c) 



T f = 




at x — Siw ? 




r)--T v 


= 


for all x £ (5, S -+ 


- 1 - f) , 


rpV 

1 wl 


= ft 
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at x = 6, 




U X J - W \ 


— d-T f 
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rpV 
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Table 3.2 

A list of dimensionless parameters and their numerical values. 



Symbol 


Definition 


Value 


Bi 


hR v 
k w 


3.45 X 10~ 4 


a 


KApi 
p w gWV v 


0.258 


P 


Pi 
Pw 


0.917 


7 
.5 


L" 
TJ 

Rf 
R-' 


0.500 
0.175 


V 


k w 


0.0419 





T a —T c 
T c 


1.83 x 10~ 5 


LU 


a 
pi?" 


0.0374 



Vessel water (right): d xx T w2 = for all x G (S + 1 + f, 6 + 2) , (3.33a) 

d x f w2 = Bi(l - Tl 2 ) at x = 6 + 2. (3.33b) 

The dimensionless parameter Bi = hR v /k w is the Biot number. 

Next are five algebraic equations for p£(i), p v g {x,t), p^,(t), c" g (t) and p g (t): 

p» g =p» g (l + 6TZ), (3.35) 
Pl=P v g -^ (3.36) 
c" g =Hp» g , (3.37) 

(f(0)) 2 -c v Jl-f 2 ) 
Pi = / (3-38) 

The non-dimensionalized Stefan condition (|3.12p becomes 

Siw — ^x^w ^ ^ = Siw 3 (3.39) 

where we have chosen the time scale 

i="% {R Zy (3-40) 

so as to eliminate the coefficient in the Stefan condition. This is clearly an appropriate time scale 
for our problem because the melting process governs the dynamics of the ice layer and hence also 
the water transport from fiber to vessel. Using parameter values from Table 13.11 we find that a 
typical value of the time scale is i ss 46.1 s. Finally, the remaining differential equations (|3.10p . 
(|3.15[) and (|3.17p for the quantities s g i(t), r(t) and U(t) become 

2/3s gi S gi + 2(1 - /3)s iw s iw - ~/U = 0, (3.41) 
2fr = -NU, (3.42) 
NU = -a [pi - p[ - (1 + Off (5, F))c v s ] ■ (3-43) 

3.5. Onset of gas dissolution in fiber after ice melts completely. As the ice in the 

fiber melts, we eventually reach a critical time when the ice layer disappears and the two interfaces 
s g i and Si W merge into a single gas-water interface whose position we denote by x = s(t). Similar 



MODEL OF SAP EXUDATION IN MAPLE TREES 



13 



to the gas/water interface in the vessel, the dynamics of s are governed by dissolution of gas within 
the fiber water and changes in pressure/volume of the gas bubble. Assuming as before that the gas 
in the fiber diffuses much faster than it dissolves, we again employ a quasi-steady approximation 
in which the dissolved gas concentration is assumed constant at any given time. Therefore, we 
replace (|3.34j) with 

pf=pf{l + 8Tf) = pf, (3.44) 
and eliminate (|3.39|) and (|3.41j) in favor of the single equation 



2ss = jU. 



(3.45) 



These last two equations are already dimensionless and are analogous to (|3.35p and (|3.42|) for 
the vessel. The modified system has as additional unknowns the water pressure pi, dissolved gas 
concentration c£, and gas density pf, which are governed by the following analogues of the vessel 
equations (I3.38|) and (|3.37|l : 



Pg- ~ 

J s 

4 = H & 

P S = =3-* = 



(3.46) 
(3.47) 

(3.48) 



where £ = Pg(0)/ pg(0) is a dimensionless parameter corresponding to the ratio of initial gas 
densities. Finally, without the ice layer to transmit the gas pressure directly to the fiber water 
located on the other side, Darcy's law (|3.43l) must be modified so as to replace the fiber gas pressure 
by pi- 



nu = -a [pi - pi - (i + erf (8, t)K] 



(3.49) 



3.6. Analytical solution for temperatures. Because the quasi-steady temperature equa- 
tions (|3.30|) - (|3.33[) have such a simple form, they can be solved analytically. In fact, the tem- 
perature is a linear function of x in each sub-domain, and applying the corresponding boundary 
and matching conditions leads to the following explicit solutions in terms of x, t and the non- 
dimensional parameters: 



T°(x,t) 

The temperature solution on the whole domain may then be summarized as 







r}(x 




■n(8 ^ 


-2 + 


Br 1 ) 4 


-2(l-t;)f(t)-tjSto(*) 


X 


-(1 


-V) (6- 


f l-f(F)) -r)s iw (i) 


T]{8 -) 


-2 + 


Br 1 )^ 


-2(l-t;)f(t)-tjSto(*) 




T)X 


+ 2(1- 




T}(6 4 


-2 + 


Bi" 1 ) 4 


- 2(1 -»j)r(t) 



(3.50) 
(3.51) 
(3.52) 



f(x,t) 



0, if x e [0, s iw ], 

ff(x,t), if x G (s iw , 8 + 1 — f], 

fv( x ,t), itxe(S + l-f,8+l + f], 

{fZ 2 (x,t), if xe (8+l + f, 8 + 2}. 



(3.53) 



Note that the matching conditions on temperature and heat flux at the vessel wall imply that the 
functions T£ and are identical. 
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3.7. Initial conditions. We choose a "base case" for our simulations in which the initial 
conditions are taken consistent with the typical scenario described by Tyree [26] : 

• Pg(0) = 100 kPa: gas in the vessel is initially at atmospheric pressure. 

• Pg(Q) = 200 kPa: gas in the fiber is at twice the pressure in the vessel owing to compression 
by the frozen sap. Based on Tyree's data, we expect that this value is an upper bound, so 
that p f g (0) actually lies in the range [100, 200] kPa. 

• s™(0) = S, s g i(0) = 6/^/2 fa 0.75 and U{0) = 0: corresponding to a fiber gas bubble 
that is initially half the total fiber volume, with the other half filled by ice. This initial 
state is established during the previous fall/winter freeze, when gas-filled fibers are at 
atmospheric pressure in equilibrium with the vessels. As the temperature falls below 
freezing, ice forms on the inside of the fiber wall and the gas pressure increases by a factor 
inversely proportional to the ice thickness according to equation ()3.34j) : hence, the value 
s g i(0) = 0.75 is consistent with the initial gas pressures chosen above. 

• r(0): is the quantity in this study with the greatest uncertainty since we haven't yet found 
data on gas concentration or bubble size in the vessels; hence, this is a natural parameter 
to vary in our simulations. We expect that the initial vessel volume fraction taken up by 
gas is significantly less than in the fiber, otherwise a large embolus is likely to remain in 
the vessel and hinder sap transport once transpiration recommences. 

Initial values are not required for the temperatures, although they can be obtained by substituting 
the above values into (|3.50[) - (|3.53|) and the algebraic relations in Section f3T4l 

4. Numerical simulations. Owing to the nonlinear coupling in the governing equations 
derived in the previous section, numerical simulations are required in order to obtain solutions 
of the full system. We employ a straightforward numerical approach in which the four ordinary 
differential equations f|3 .39[) and (|3.41[) - (13.43[) are solved for s g i, S{ W , f and U using the stiff solver 
odel5s in MATLAB, with error tolerances RelTol = 10~ 10 and AbsTol = 10~ 8 . Temperatures 
are calculated using the analytical expressions in (|3.50[) - p.52p . and the system is closed using the 
algebraic equations p.34p - (|3.38[) . 

The above procedure applies when ice is present in the fiber. In order to permit MATLAB to 
modify the equations automatically once the ice disappears, the event detection feature in ode 15s 
is used to determine the instant when s g i = S{ W and V/ — > 0; at this time, the simulation is 
halted and the calculation is restarted with the equations modified as described in Section [331 To 
assist in understanding which equations are solved under which circumstances, we have provided 
in Figure |4~T1 a graphical summary of the equations and variables in the two cases. 

In the remainder of this section, we perform a series of simulations with increasing model 
complexity, first studying the dynamics of gas dissolution and flow through the fiber wall when 
the pressurized fibers contain no ice, then introducing osmotic effects, and finally incorporating 
the ice layer and phase change. Even though the code is written in dimensionless variables, all 
parameters and results from this point on are converted to dimensional variables using equations 

(EUD. 

4.1. Case 1: Dissolution and porous transport only (no ice, no osmosis). We begin 
by considering the simple situation in which there is no sucrose in the vessel and no ice in the fiber 
(and hence, no phase change). As a result, the sap dynamics are governed primarily by two effects: 
gas dissolution and the flow of sap driven by pressure exchange between gas bubbles in the fiber 
and vessel. We note that temperature variations are relatively small and so they have a negligible 
effect on the solution. 

Gas bubble dissolution has been studied extensively in the context of sap flow by several au- 
thors [lOl EU [33] , who concluded based on experimental evidence and simple physical arguments 
that pressurized gas bubbles in xylem sap should dissolve completely over a time period ranging 
anywhere from 20 minutes to several hours. These results are partially supported by the mathe- 
matical study of bubble dissolution by Keller [9] , who proved that for an under-saturated solution 
(such as xylem sap) that fills an unbounded domain, the only stable equilibrium solution is one in 
which the gas bubble shrinks and eventually disappears. An analogous situation was studied by 
Yang and Tyree [33], who performed experiments on excised maple branches open to the air and 
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Fig. 4.1. List of variables and governing equations depending on whether or not ice is present in the fiber. 



developed a corresponding mathematical model in which the pressure is maintained at atmospheric 
conditions at the open boundary. Their experimental and numerical results are consistent with 
the complete dissolution of gas bubbles in finite time. 

In our study of sap exudation however, we are interested in the behavior of xylem sap within 
an intact maple tree stem that is not exposed to the open air; consequently, the open boundary 
conditions treated by Yang-Tyree above are not applicable here. Fortunately, Keller also considered 
a pressurized gas bubble suspended in a closed fluid container [5] for which he proved that the 
bubble will still shrink in size as it dissolves, but in this case will reach an equilibrium state having 
non-zero radius. This behavior is contrary to that observed above in experiments on excised maple 
branches, and so bears careful consideration. 

To investigate the ability of our model to capture Keller's results for a closed container, we 
performed several simulations with bubbles of various sizes in both fiber and vessel. We start 
with a scenario in which the fiber bubble takes up 50% of the fiber volume and the vessel bubble 
only 10%, corresponding to s gl (0) = R f /V2 « 0.7 Rf and r(0) = R v /VW « 0.3R V . Based on 
the discussion of initial conditions in Section 13.71 the gas in the vessel is initially at atmospheric 
pressure (100 kPa) while that in the fiber is twice as large (200 kPa). The resulting dynamics 
for the pressures and bubble radii are shown in Figure 14.21 from which we observe that the fiber 
bubble grows in size while the vessel bubble shrinks, until both reach an equilibrium state at a time 
of roughly t = 40 seconds. Figure B~27 a) shows shows that the system approaches an equilibrium 
state in which the water pressure in the two compartments is equal (at approximately 155 kPa). 
This is clearly consistent with (|3.43|) since there can be no flux across the porous cell wall at the 
equilibrium state, which requires that p{, = p-^ (remembering also that 5^=0). The plots show 
that the gas and water pressures vary inversely with the radii as expected, and the piecewise linear 
dependence of temperature on position in equations (|3.50l) - (|3.53[) is evident in Figure H^d). 

The plot of gas pressure shows a decrease in pf over time as the fiber bubble expands and 
a corresponding decrease in p v as the vessel bubble is compressed. This behavior is consistent 
with Keller's results for bubbles in a closed system, and also captures the transfer of pressure from 
fiber to vessel that we would expect in sap exudation. It is important to note that the increase 
in vessel water pressure from roughly 90 to 155 kPa (see Figure B~2l a)) is close to the upper limit 
of the range of 30-60 kPa that is to be expected for positive pressures in maple according to [25] . 
Finally, we note that the magnitude of the pressure change in the fiber arising from dissolution of 
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Fig. 4.2. Case 1 - base case: No ice, initial radii s(0) = 0.7Rf , r(0) = 0.3R V , and pressure p g {0) = 200 kPa. 



gas is only 5 kPa, which is small relative to the other pressure differences between fiber and vessel; 
consequently, dissolution has only a small effect on sap flow over this time scale. 

The parameter value with the most uncertainty is the initial value of vessel bubble radius, 
and so we focus next on the effect of reducing r(0) in Figure l4~3l As long as r(0) is greater than 
approximately 0.2R V , we observe behavior similar to the earlier simulations in Figure l4~2l in that 
the fiber gas bubble expands, driving water through the cell wall and hence compressing the gas 
in the vessel. However, when r(0) = 0.2R V the vessel gas bubble dissolves entirely, causing the gas 
pressure to become unbounded as r — > owing to the air/water surface tension term in (|3.36[) . 
At the precise moment when the vessel bubble disappears, the model breaks down because there 
is no further transfer of pressure between fiber and vessel, nor is there any flow of sap through 
the fiber/vessel wall. This point actually corresponds to a terminal steady state of the system, 
and in our computations we have found it necessary to specify a very small threshold value of the 
bubble radius below which we halt the simulation. The collapsing vessel bubble case r(0) = 0.2R" 
is illustrated in more detail in the plots in Figure 14.41 

The dissolution of gas in the vessel has an intriguing connection with the phenomenon of 
embolism repair |33j . in which gas bubbles within the vessel that could potentially hinder sap 
transport are known to dissolve during the spring thaw. The physical mechanism underlying 
embolism repair is still not fully understood, and has implications for a wide range of other tree 
species; hence, a more detailed study of bubble dissolution as relates to embolism would form a 
fascinating avenue for future research. 

Since pf (0) = 200 kPa is at the upper end of the range of fiber gas pressures we would expect 
in actual trees, we next consider the effect of varying p^(0) over the range [100,200] kPa. As 
seen in Figure 14.51 the qualitative behavior is similar until the initial fiber gas pressure approaches 
the atmospheric level in the vessel; for pl(0) low enough, the dynamics are reversed in the sense 
that the fiber bubble shrinks and the vessel bubble expands. As seen in the detailed plots for 
pl(0) = 100 kPa in Figure l4~6l the fiber and vessel again equilibrate at a state where the water 
pressures are equal, although in this case the level of p w sa 75 kPa is much lower than that 
observed in the base case. Clearly, when the fiber begin at lower pressure, then the tree is less able 
to pressurize the vessel sap; furthermore, the fiber compartment must be above some threshold 
pressure for there to be any pressure transfer from fiber to vessel that is consistent with exudation. 

In fact, if the initial fiber gas pressure and radius are both small enough, then it is the gas 
in the fiber that dissolves entirely. Although this situation is not likely to occur in practice, we 
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Fig. 4.3. Case 1 - varying r(0); No ice, initial radius s(0) = 0.7 R$ , and pressure Pg{0) = 200 kPa. 
(a) Water pressures (b) Gas pressures (c) Gas bubble radii 
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FlG. 4.4. Case 1 - collapsing vessel bubble: No ice, initial radii s(0) = 0.7 R$ , r(0) = Q.2R V , and pressure 
Pg(0) = 200 kPa. Here, the initial vessel bubble is so small that it dissolves entirely. 
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Fig. 4.5. Case 1 - varying p g (0): No ice, initial radii s(0) = 0.7 R f , r(0) = 0.3R". 
(a) Water pressures (b) Gas pressures (c) Gas bubble radii 
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FlG. 4.6. Case 1 - reduced fiber pressure: No ice, initial radii s(0) = 0.7W and r(0) = 0.3H", and pressure 
p f g (0) = 100 kPa. 
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have nonetheless included a simulation with p£(0) = 100 kPa and s(0) = 0.25i?^ in Figure I4T41 to 
illustrate this point. 
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Fig. 4.7. Case 1 - collapsing fiber bubble: No ice, initial radii s(0) = 0.25-R-f , r(0) = 0AR v , and pressure 
Pg(0) = 100 kPa. The fiber bubble radius and pressure are initially so small in this case that the fiber bubble 
dissolves entirely. 

4.2. Case 2: No ice, with osmosis. For the next round of simulations, we introduce 
osmotic effects arising from the presence of sucrose in the vessel, and still maintain our focus 
on the effect of porous flow and gas dissolution by assuming all sap is in the liquid phase. The 
vessel sap is taken to contain 2% sucrose by mass (with 0% in the fiber), but otherwise the initial 
conditions are similar to Case 1. 

The first simulation with r(0) = 0.3R V , p v g (0) = 100 kPa and p f g (0) = 200 kPa is displayed 
in Figure S^H] and should be compared with Figure |4~21 (base case, without osmosis). Osmotic 
pressure has a relatively small impact on the fiber bubble dynamics, with the fiber radius being 
only slightly smaller and the equilibrium water pressure dropping from 155 kPa to just under 
150 kPa. In contrast, there is a major change in the vessel bubble evolution, with the equilibrium 
radius dropping to roughly half that of the base case. The source of the difference is evident from 
the pressure plots: the initial vessel water pressure is nearly 140 kPa lower than the base case in 
Figure U^a) owing to the osmotic pressure, and the equilibrium vessel gas pressure increases by 
roughly the same amount. 

Maintaining the fiber bubble at pressure Pg(0) = 100 kPa, we next draw comparisons with 
Figure l4~3l which varied the vessel bubble radius over values r(0) = {0.2,0.25,0.3,0.4}^. From 
the summary of results shown in Figure I4T51 we see that the presence of osmosis enhances the flow 
from fiber to vessel and leads to a much larger pressure increase in the vessel sap. The smallest 
vessel gas bubble with r(0) = 0.2 disappears sooner in the presence of osmosis, which is to be 
expected. 

A more insightful comparison is afforded by reducing the fiber gas pressure to Pg(0) — 100 kPa, 
and comparing Figure 14.101 to the corresponding case without osmosis in Figure 14.61 (where we 
observed a reversed sap flow, from vessel to fiber). Clearly, osmotic effects act to return the sap 
dynamics to "normal," with water now being driven from fiber to vessel. Except for the shift 
in equilibrium pressures, the qualitative solution behavior is now consistent with the base case 
simulation in Figure 14.21 Consequently, our model demonstrates that osmotic pressures owing to 
sucrose in the vessel can provide a mechanism for enhancing sap exudation, even when the gas 
pressure in the fiber is as low as that in the vessel. In this respect, our model is consistent with 
the hypothesis of Tyree [2 6) that claims osmosis prevents gas bubble dissolution in the fiber and 
sustains sap exudation pressures over longer time periods. 

4.3. Case 3: Milburn O'Malley model (with ice, no osmosis). We next test Milburn 
and O'Malley's hypothesis for sap exudation that includes the effects of melting ice and gas disso- 
lution, but ignores osmosis. We restrict ourselves to the time period over which ice is present, so 
that gas in the fiber is isolated from the water phase by the frozen ice layer and hence gas disso- 
lution occurs only in the vessel (the long-time dynamics after the ice layer melts are considered in 
Case 5). Furthermore, gas-water interfacial tension plays no role in the fiber in this case, so that 
the pressure of the melted liquid adjacent to the fiber wall is simply equal to the gas pressure. 
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Fig. 4.8. Case 2 - base case with osmosis: No ice, with osmosis, initial radii s(0) = 0.7 W , r(0) = 0.3R? , and 
pressure Pg(0) = 200 kPa. 
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Fig. 4.9. Case 2 - varying r(0): No ice, with osmosis, initial radius s(0) = 0.7 , and pressure Pg(0) = 200 kPa. 



In Figure l4~.lH we present the simulation results using initial conditions s g i = 0.7R*, s W i = Rf, 
r(0) = 0.3R V , pjj(0) = 200 kPa and Pg(0) — 100 kPa, which can be compared with the base case 
depicted in Figure B~2"1 Until a time of roughly 40 s, the qualitative features of the two solutions 
are similar, with the water pressure decreasing in the fiber and increasing in the vessel until they 
equilibrate at some intermediate value. Over this time period, the fiber bubble radius increases 
smoothly from the initial 2.45 /zm to roughly 2.6 /im, and the motion of the vessel bubble (not 
shown) is also similar to that of the base case. There is, however, a significant difference in the 
gas and water pressures owing to the absence of surface tension. 

After the initial rapid equilibration phase, the two water pressures remain equal and experience 
a gradual decrease over a much longer time scale. This slow adjustment in pressure is due to 
phase change and is attributed to the roughly 10% difference in density between water and ice. 
The ice/water interface Si W (t) has a roughly linear variation in time according to Figure FLllf c). 
which should be contrasted with the \[i dependence that typically arises in Stefan problems; this 
discrepancy is due to the fact that sap outflow through the fiber wall dominates the interface 
motion. The intersection of the s g i and Si W curves near t ks 2000 s represents the time when the 
ice layer disappears. Our main conclusion from these results is that even when ice is present, the 
vessel compartment is pressurized by the fiber, with the main difference (comparing Figures 14.21 
and 14. lip being that the vessel sap pressure is higher because of the missing interfacial tension 
force. 

A final series of simulations is then performed in which the initial vessel bubble radius is varied 
over the range r(0) = {0.2, 0.25, 0.3, 0. 4}R V . Figure H. 121 shows that the smallest bubble dissolves 
entirely within approximately 10 seconds, after which the transfer of pressure from fiber to vessel 
stops. In all other cases, the fiber and vessel compartments equilibrate over a time scale of roughly 
2 hours. 
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Fig. 4.10. Case 2 - reduced fiber pressure: No ice, with osmosis, initial radii s(0) = 0.7 R? , r(0) = 0.3R? , and 
pressure Pg(0) = 100 kPa. 
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Fig. 4.11. Case 3 - base case with ice: Milburn-O 'Malley scenario, no osmosis. Initial radii are s g i(0) 
0.7R f , s iw = R f , r(0) = 0.3R V , and pressure p f g (0) = 200 kPa. 



4.4. Case 4: Full model, with ice and osmosis. In this section, we simulate the full 
model including ice and osmosis, beginning again with the base case. The solution plots displayed 
in Figure 14.131 should be compared with the Milburn-0 'Malley simulation in Figure [4.111 which 
differs only in the lack of osmosis. The two sets of results are qualitatively similar, there being 
slight changes in the water pressure and bubble radii owing to the inclusion of osmotic pressure. 
Consistent with the ice-free results in Cases 1 and 2, the vessel gas pressure is increased by an 
amount roughly equal to the osmotic pressure II = c s lZT « 132 kPa. 

If the initial vessel bubble radius is decreased slightly from 0.30H" to 0.29R V , then the vessel 
bubble is no longer unable to equilibrate with the fiber and dissolves completely as shown in 
Figure 14.141 This is consistent with the other results in Figures 14.41 and 14. 9[ except that the extra 
osmotic pressure in the vessel causes the gas to disappear at larger values of r(0). When viewed 
in the context of embolism recovery, this is evidence that suggests osmosis can enhance the ability 
of the tree to recover from embolism. 

A series of simulations for different values of the initial vessel bubble radius are reported in 
Figure |4.15[ which shows a similar outcome as in Case 3. The main difference is that the bubbles 
with radii r(0) = 0.2i?" and 0.25i? u both dissolve completely, which is further evidence that osmosis 
enhances the effect of exudation and makes the vessel bubbles more likely to collapse. 

4.5. Case 5: Long-time behavior after ice melts completely. The purpose of this final 
set of simulations is to study the behavior over the longer term once the ice is totally melted, and 
to demonstrate that our approach is capable of computing through the disappearance of the ice 
layer and the corresponding change in the governing equations. We present results for a few choices 
of parameters considered in the previous section (Case 4) and emphasize that until the time that 
the ice melts, the plots are identical. Our main observation is that after the ice disappears, the 
solution dynamics exhibit a rapid equilibration to a constant equilibrium state that is consistent 
with what we saw earlier in Case 2 (no ice, with osmosis). 

Simulations are shown in Figures [4.161 and 14 . 1 71 for two values of vessel radius, r(0) = 0.3 and 
0.6, that can be compared with results in Figure 14.131 and f4. 151 At the instant the ice layer melts 



MODEL OF SAP EXUDATION IN MAPLE TREES 
(a) Vessel bubble radius (b) Vessel water pressure 



21 















200 
















5 10 


10 
















10 

9 








180 














S 






5 










7 










160 





5 10 


160 






6 



































140 
120 


-- r =0 - 2RV 




5 
4 




Pressure 


140 






i 


r =0.25R v 




3 






120 






100 


--- r o =03RV 




2 
1 




r^O.ZR" r =0.25R v 


r Q =0.3R v r =0.4R v 




100 










r Q =0.4R v _ 






( 


























1000 2000 3000 4000 5000 
time (seconds) 









1000 2000 3000 
time (seconds) 


4000 5000 



Fig. 4.12. Case 3 - varying r(0): Milburn-O'M alley scenario with ice and no osmosis, initial radius s(0) 
0.7 Rl , and pressure Pg(0) = 200 kPa. The "zoomed-in" sub-plot shows the curves near t = 0. 



(a) Water pressures 



(b) Gas pressures 



(c) Fiber phase interfaces 




Fig. 4.13. Case 4 - base case with ice and osmosis: Initial radii s 9 ;(0) = 0.7W , r(0) = Q.3R V , and pressure 
p f g (0) = 200 kPa. 



away, the water pressures undergo a sudden transition in which and p^, drop sharply and then 
rapidly equilibrate to a new constant value. This pressure discontinuity arises from the sudden 
appearance of the surface tension term in equation (|3.46p , when the ice layer vanishes and the fiber 
sap comes into contact with the gas. We emphasize here that the numerical algorithm handles 
the discontinuity and change in governing equations easily and without introducing any spurious 
oscillations or other non-physical effects. 

After the jump in water pressure, there is a rapid but smooth variation in the radii and gas 
pressure as the system equilibrates (although in the case when r(0) = 0.3R V the change in gas 
pressure is too small to be seen on this scale). Otherwise, the behavior after the loss of ice is 
analogous to the earlier ice-free simulations. 

In reality, we would expect to see a rapid but smooth variation in liquid surface tension as the 
disappearing ice layer transitions through a "mushy region" where both ice and water coexist [?]. 
Including this additional level of detail would not have a major effect on our solution, but it could 
form the basis for an interesting future study of the detailed structure of the transition layer. 

5. Conclusions and future work. The aim of this work was to develop a detailed mathe- 
matical model that captures the physical effects - heat and mass transport, gas dissolution, surface 
tension, osmosis and phase change - that are believed to play a significant role in sap exudation 
in maple trees. In particular, we were interested in investigating two competing hypotheses for 
sap exudation based on a model of Milburn-O'Malley (without osmosis) [12] and another of Tyree 
(with osmosis) [26] . The model focuses on the cellular scale and the transfer of pressure between 
two main xylem components: fibers and vessels. We derived a coupled system of nonlinear, time- 
dependent, differential-algebraic equations, which we solved using MATLAB. Extensive numerical 
simulations were performed that uncovered the following: 

• It is necessary to include the effect of gas bubbles in the vessel sap in order to allow a 
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Fig. 4.14. Case 4 - collapsing vessel bubble: With ice and osmosis, initial radii s gt (0) = 0.7 R f , r(0) = 0.29K", 
and pressure p g (0) = 200 kPa. 
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Fig. 4.15. Case 4 - varying r(0): With ice and osmosis, initial fiber radius s g i(0) = , and pressure 

Pg(0) = 200 kPa. The endpoint of each curve corresponds to the time that the ice layer completely melts. The 
"zoomed-in" sub-plot shows the curves near t = 0. 
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(b) Gas pressures 



fiber 

vessel 



1000 2000 3000 4000 5000 
time (seconds) 



= 4 

i 3 

- 2 
1 



(c) Gas bubble radii 



1000 2000 3000 4000 5000 
time (seconds) 



Fig. 4.16. Case 5 - base case with ice and osmosis: Long-time behavior with initial radii s g i(0) = 0.7R^ , 
r(0) = 0.3R V , and pressure p f g (0) = 200 kPa. 
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Fig. 4.17. Case 5 - larger vessel bubble: With ice and osmosis, initial radii s gl (0) = 0.7 R f , r(0) = 0.6R V , 
and pressure p g (0) = 200 kPa. 
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transfer of stored pressure in the fibers to the vessel sap. This gas in the vessels was not 
explicitly mentioned by either Milburn-O'Malley or Tyree. 

• Taking parameter values consistent with data in the literature, our model captures positive 
vessel pressures provided that the gas in the fiber is sufficiently compressed during the 
previous year's freezing phase. The observed increase in vessel sap pressure lies in the 
range of 30-60 kPa that is consistent with experimental values reported in the literature 
for maple sap exudation. 

• For a wide range of parameters, and provided the vessel gas content is low enough, the fibers 
are able to completely dissolve the gas bubbles within the vessel sap. This is consistent 
with the phenomenon of embolism repair that is known to occur in maple as well as other 
hardwood species. 

• Reasonable exudation pressures are obtained both with and without osmosis. One impact 
of including osmotic effects is to increase slightly the threshold vessel bubble size below 
which gas in the vessel gas dissolves entirely; in this respect, our model suggests that 
osmosis may enhance embolism repair. The second main impact of osmosis is to lessen the 
likelihood that the fiber gas bubble collapses, which enhances the ability to exude sap. 

This paper can only be considered as a preliminary study of sap exudation because the com- 
parisons and conclusions we have drawn so far are mainly qualitative. More detailed comparisons 
with experimental data in the literature are currently underway to further validate the model, and 
will form the basis of a future publication. There are a number of other specific areas that we 
believe would provide fruitful avenues for future research: 

• The long-time behavior pictured in Figures 14.161 and 14.171 suggests that this problem is 
ripe for a multi-layer asymptotic analysis in which the solution is divided into two pieces: 
a linear solution, matched to initial conditions via a boundary layer on the left; and a 
constant equilibrium state that is matched across an interior layer to the linear solution. 

• The predictive power of this model would be strengthened with better estimates of two 
parameters: the gas content in the vessel, and the initial pressure in the fiber. We will 
perform a more extensive search of the literature and focus in particular on data from 
studies of embolism. In conjunction with this effort, we believe that our model (with only 
minor modifications) could also be applied to study the phenomenon of embolism repair 
coinciding with the spring thaw 33 . 

• The conditions under which the maple tree freezes at the end of the previous season are 
known to have a significant impact on sap exudation during the following spring. To study 
this process, we could run our model "in reverse" to study the cooling sequence in the 
Milburn-O'Malley model (steps 1-4 in Figure [2T2|) . 

• We are currently developing a macroscopic tree-level model, based on the ID model of 
Chuang et al. [3] that uses Richards' equation to capture sap flow in the porous xylem. 
In order to apply this work to the study of exudation, we require solution-dependent 
transport coefficients that will be derived from the microscopic cell-level model derived in 
the current paper using an up-scaling procedure. 
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